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ABSTRACT 

A recursive algorithm was designed to invert large, dense, symmetric, positive- 
definite matrices using small amounts of computer core, i.e, , a small fraction 
of the core needed to store the complete matrix. The algorithm described in 
ttiis report is a genei'alised Gaussian elimination technique. Other algorithms 
are also discussed for the Cholesky decomposition and step-inversion techniques. 

The xairpobe of the inversion algorithm is to solve large linear systems of nor- 
mal equations generated by worldng geodetic problems. The algorithm was in- 
corporated into a computer program called SOLVE (Reference 1 and 2). In the 
past the SOLVE program has been used in obtaining solutions published as the 
Goddard Earth Models (References 3 through 6). The latest of the Goddard 
Earth Model (GEM) solutions, GEM S, contained approximately 1156 modeled 
variables. 
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SECTION 1 - INTRODUCTION 

This report is concerned with the solution of linear systems of least-squares 
normal equations. The problem, described in matrix notation as 

Ax = y 

is to solve the vector, x , given the vector y and a normal matrix, A , 
(positive-definite and symmetric). 

The matrix inversion algoritlim described in this report was designed for in- 
verting large matrices using a fraction of the computer storage space (core) 
needed to store the complete matrix. The algorithm was designed to partition 
the matrices into small segments which fit into the available computer core. 

Section 2 describes a generalized form of Gaussian elimination. This algorithm 
is presently used in the SOLVE 2 (Reference 1) program at Goddard Space Flight 
Center (GSFC) and is being programmed into several other geodetic computer 
programs as well. The algorithm is programmed in a package which includes 
a subroutine, PINV, which is the main driver. The PINV package is presented 
in the Appendix. 

Section 3 describes two alternative algorithms using the same basic partitioning 
technique described in Section 1. These algorithms are Cholesky decomposition 
and a step-inversion method. 

Section 4 deals with numerical problems in obtaining a matrix inverse and pre- 
sents a simple metliod of obtaining a set of "condition niunbers" which indicate 
numerical loss of precision for each row of the inverse matrix. Also described 
is a simple algorithm to iterate a solution vector in order to minimize numerical 
truncation error. 

Test results using the PINV inversion paclcage are presented in Section 5. 
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SECTION 2 - A GENEHiVLIZED GAUSSIAN 
EIJMINATION ALGORITHM 


This section presents an algorithm for inverting symmetric, positive-definite 
matrices using a generalized form of Gaussian elimination. Normally in the 
elimination process, one row is eliminated at a time. The generalized algo- 
rithm allows for a set of rows, called a segment, to be eliminated each time. 
The size of the segment is selected to fit into available computer core, v/Mle 
the remainder of the matrix is stored on external storage devices. The algo- 
rithm. presented in this section is coded in FORTRAN and listed m the Appendix. 

The equation 


Bx - y 


can readily be partitioned to take the form 



From the symmetry properties of B 


T 


22 22 


®21 = 


OBlGIN’Aij PACF to 
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(2-1) 


(2-2) 

(2-3) 

(2-4) 

(2-5) 





The size of B can be selected so that B and B contain no more ele- 
XX XX x^ 

ments than the computer storage will allow. This usually means that B will 
not fit into the same storage area and, therefore, will need to be recursively 
partitioned in the same manner as B . 

If B is nonsingular, a matrix, V , exists such that 


^11 

^21 ^ 22 / V 2 / 


(2-6) 


and 


'=11 =12W''ll ^12\ 

'=21 = 22 ) ki O 'o I 


< 2 - 7 ) 


In order to compute the submatrices of V , intermediate matrix sets are de- 
fined by 


C = B ^ 
11 11 


(2-8) 


^12 = ^11 ®12 


( 2 - 9 ) 


^21 


(2-10) 


22 22 " 21 12 


( 2 - 11 ) 
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The solutions for the submatrices of V are then 


V = C 
22 22 


V = -C V 
12 12 22 


T 

V = V 
21 12 


V = C - V c 
11 11 12 21 


Solutions for the vector x are 


(2-12) 

(2-13) 

(2-14) 

(2-15) 


*2 = °22 [^2 - ‘'21 ^l] 


*1 ^11 ^1 “ ^12 *2 


(2-16) 

(2-17) 


A more direct method of computing the solution vector x would be to use the 
entire inverse matrix, i. e, , 

x=Vy (2-18) 

However, the advantage of using the generalized Gaussian method (Equa- 
tions (2-16) and (2-17)) is that the solution can be obtained without the complete 
inverse matrix. 

Given below is a recursive algorithm for this method in which the matrix B 
is stored on an external scratch device and submatrices read into the internal 


core storage area as needed. Only the upper triangular portion of B elenaents 
(on and above the diagonals) need be stored because of the symmetric property 
of B . Three scratch units, denoted as units 1, 2, and 3, are used in the algo- 
rithm. The steps in the algorithm are as follows: 

1. Read as many rows of B on xmit 1 as possible to fill the available 
core storage area. The submatrices B, - and B are formed 

11 1a 

from these rows. B is also initiated into core by using the sym- 

^1 rjp 

metry property B„ = B (see step 2 below). 

a1 1a 

2. Compute 



‘^12 = ^12 


^2 = ^ 2 ' °21 ^1 


and 


v' = C V 
•^1 11 


simultaneously replacing B » B , y , y in core with the 

11 1a 1 a 

newly computed variables, C.. , C _ , y' , y’ . 

11 1a 1 a 

3. Store the and submatrices on a scratch storage device 

(i. e. , unit 3), 

4. Read an additional row of the matrix B from scratch unit 1 which 
belongs to the submatrix B on unit 1. Compute one row of the 
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C submatrix (Equation (2-11)) and write the row on unit 2. For 
the ith row of , the equation for the jth element is 


C22(i.j) = ^ ‘ 

k 


Repeat step 4 until all rows of the B matrix have been processed. 

Rewind scratch units 1 and 2. The identification numbers for the 
two scratch units are switched so that the C matrix can be re- 
cursively partitioned in the same manner (i. e. , steps 1 through 4 
are repeated) until the matrix written on scratch unit 2 is small 
enough to fit into the available core area. (Fig'ure 2-1 illustrates 
the information on the external scratch units for steps 1 through 5. ) 


Read from unit 2 and store in core. Compute 

and write the results back on unit 1. Rewind unit 2. Compute 


*2 “ '^22 5'2 • 


Read the last C . , C submatrix sets written on unit 3 and store 
XI 1 a 

in core. Read through unit 1; operating with one row of at a 

time, compute = “^12 ^22 * 


Compute , where 


Si 


and add the V_ , and V. „ submatrices with the V matrix on 

11 1a Za 

unit 1. Rewind unit 1. If the complete matrix B has not been in- 
verted, return to step 7. (Figure 2-2 illustrates the data stored 
on the scratch units during steps 6 through 8. ) 
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UNIT 1 


UNIT2 


UNIT3 



V 1 

^12 

N 



Figure 2-1. Information Stored on Scratch Units in 
Steps 1 Through 5 for the Recursive 
Gaussian Method 


UNIT1 UN1T3 



Figure 2-2. Information Stored on Scratch Units in 


Steps 6 Through 8 for the Recursive 
Gaussian Method 
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SECTION 3 - CONSIDER/iTION OF OTHER INVERSION METHODS 


The basic idea of segmenting the rows of the matrix to be inverted can be ap- 
plied to other matrix inversion methods. This section describes two alternate 
inversion algorithms. 

3. 1 CHOLESKY DECOMPOSITION BY PARTITIONING 

Equation (3-1) can be reformed by assuming the existence of a matrix, C , such 
that 


C^C = B (3-1) 


and 


T 

C Cx=y (3-2) 

where the matrix C has elements whose values are zero below the diagonal. 
Matrix C is the decomposition of B . The vector w is then defined such that 

w = Cx (3-3) 

and, therefore, 

c'^w = y (3-4) 

T 

Since the elements of C are zero above the diagonal, w can be solved using 
forward substitution. Once w is obtained, x can be solved using back sub- 
stitution. 
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Since C has values of zero below the diagonal, 


c. . = b.. 

11 1 ] ]1 


(3-5) 


B is symmetric; therefore, 


b„ = b . 
]1 1] 


(3-6) 


“ii N 


(3-7) 


and 


“ii “ \i 


(3-8) 


Therefore, in order to decompose the first row of B , the elements of C are 




ij yb. 


11 


(3-9) 


The matrices C and B can be partitioned so that 


° \Pll % ®12 


<2 ‘^ 22/\0 ^ 2 , 



(3-10) 


®21 ® 22 / 
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From, this equation, 


^12 ^12 ^22 ^22 ®22 


(3-11) 


By partitioning so that C and C are the first row of the matrix C , and 

1.x JL/d 

by rearranging Equation (3-11) so that 


^2 S2 ^22 " ^12 ^12 


(3-12) 



can be computed b5r decomposing the new ma,trix B' , where 





(3-13) 


Using the above recursive technique, C can be computed one row at a time. 
Each time a row of C is computed, it can simply replace the same row of B 
in the computer storage. If not enough storage is available for the complete 
matrix, then only those rows which will fit into the area need be operated on. 
Once those rows are decomxDosed, Equation (3-13) is used to obtain a reduced 
matrix .B' , which can then be decomposed. 

An inverse ma,trix can be obtained by setting y in Equation (3-2) to be a col- 
umn of the identity matrix and then solving for one row at a time. The Cholesky 
method is best used for obtaining solutions v/hich do not require inverse ma- 
trices, 

3.2 RECUKSIVE STEP METHOD 

The recursive step method uses the same formulas as the generalized Gaussian 
inversion process but differs in the manner in which the formulas are employed. 
This method is designed to work more effectively using’ the lower triangular 


ORIGINAL’ PAGE m 
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matrix, whereas the Gaussian elimination algorithm uses the upper triangle. 
The inverse of one submatrix is generated on a scratch unit and the other sub- 
matrices stored in core to be used to complete the inverse. The procedure is 
recursive in that the complete inverse of a smaller matrix is the inverse of the 
submatrix of a larger matrix. The procedure is repeated in order to obtain 
the complete inverse of the larger matrix. The steps in the algorithm are as 
follows: 

1. Read from a scratch unit 1 as many rows of the matrix B that can 
be stored into the available machine core. This unit will then con- 
tain the lower triangular elements of the matrix B , designated 

submatrix B,. . 

11 

2. Invert the submatrix and write the resultant inverse 

fc,, =B,M on scratch unit 2. Rewind unit 2. 

\ 11 11 / 

3. Read from scratch unit 1 as many additional rows as possible and 
form the submatrices B and B from these rows. 

ZX did* 

4. Read one row at a time of from unit B and compute 


inserting the C matrix in the computer core. Rewind unit 2. 


5. Compute 


^22 ®22 ” ®21 *^12 
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These submatrices can be computed completely in the computer 

core area by replacing the area used for i>y and finally 

by V ; the B submatrix is then replaced by V . 

22 2X 2X 

Read through unit 2 again one row at a time; compute 


- °12 


and write the results on scratch unit 3. Rewind unit 2. 

Write the submatrices V , V on unit 3 so that the complete 

2X 22 

inverse is on unit 3. Rewind unit 3. If the inverse of the complete 
matrix B has not been computed, switch units 2 and 3 so that 
is now the matrix V ; return to step 3. (Figure 3-1 illustrates 
data on the scratch units as they are used in the above steps.) 

Compute the solution vector by reading the inverse matrix V on 
unit 3. Solve x = Vy . 



SECTION 4 - ERROR ANALYSIS 

The source of error considered in this section is the numerical round-off error 
which occurs in floating point arithmetic when obtaining the solution to the linear 
system 


Bx= y 


(4-1) 


If the assumption is made that the inverse matrix V can be obtained with the 
same numerical accuracy as the original matrix B , then the round-off error 
would occur in the matrix product 


X = Vy 


(4-2) 


The numerical round-off error should occur primarily in the process of sub- 
traction; i. e. , if two six-digit numbers, a and b , are subtracted and yield 
a four-digit difference, then the answer would be good to two less decimal 
places. For a series summation, such as 


X. = y. 

1 V 1] 3 


(4-3) 


a condition number, , is defined as 


E 




‘^1 E (X.) 


max 


(4-4) 


where E (v.. y.) is the expected value of the largest term in the summa- 
' 1 ] *^ 3 'max 

tion, and E (x^) is the expected value of the solution^ The logarithm of the 


4-1 



condition number should then yield the number of decimal places lost in the 
matrix product. 


There is a round-off error, however, that is accumulated when the V matrix 
is computed. Therefore, the condition number noted above indicates the best 
possible answer, i.e. , when V is accurate to as many decimal places as the 
original matrix B . 


In order to obtain the condition number, Cj^, the expected value of a random 
element of a vector is defined as the square root of the variance, i.e. , 

n “1 1/2 


E (X.) = 

1 


H X. 
i = l 3 


n.. 


(4-5) 


The above definition is useful only if it is reasonable to expect any one element 
of a vector to be the same magnitude" as another. We can define such a vector 
as being 'hormaly distributed”. 

Given a vector X and a normaly distribu^^?d vector y the maximum expectation 
of the vector product, as used in equation (4-4), would be 


E (X . Y, )_„ = E (X. ) ' E (Y. ) (4-6) 

' i i 'max ' i 'max ' 1 ' ' 

The above equation cannot be applied directly to equation 4-4 because the vector 
is not, in general least squares problems, normaly distributed. 

The expected magnitude of the sum of elements of a vector can be defined as 

n 

E ( i; X,) 
i»l 
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This definition is consistant with statistical error estimates of a sum of obser- 


vations. 

Equation 4-7 can then be used to define the expected magnitude of a vector prod- 
uct given the vector x and a normaly distributed vector y 

n n 

E ( E X. Y. ) = E ( X; X.) • E (y.) (4-8) 

i=l i=l 

In order to evaluate the expected values in Equation (4-4), some insight into 
the least-squares process is needed. The matrix B is obtained from a set of 
observation equations in the matrix notation 

Ax = r (4-9) 

In the equation noted above, the vector r contains a set of weighted residuals 
whose expected values are 


E (r.) - 1 , 


i.e. , CT^= ±1 for weighted observations and therefore r is considered normaly 


distributed. In the least-squares process, 


T 

B = A A 


(4-10) 


and 


aT 

y = A r 


(4-11) 


From Equation (4-11), the ith element of the vector y is 


y. = / a., r. 
•"i ij J 


(4-12) 
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and using definition 4-8 the expected value of y. is 


E (Yj) = E a. . E (r.) 


using definition 4-7 it is evaluated as 


E 




1/2 


From Equation (4-10), the right hand term is 




1/2 


and Equation (4-14) becomes 


E (y.) = (b..) 

1 u 


1/2 


From Equation (4-16), 


E (y.) E (y.) 


' 11 ]] 


for all values of i and j . If the matrix D is defined such that 


11 


(4-13) 


(4-14) 


(4-15) 


(4-16) 


(4-17) 


(4-18) 
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then Equation (4-1) would become 


D 


Xj (bj^) 


1/2 


X (b ) 
n nn 


1/2 



(4-19) 


From the positive definite properties of the matrix B and using Equation (6-18) , 


d.. 1 <1 when i / j 
ij 


= 1 when i = j 


The maximum expected value of the ith row of D , E ( d..| , is the 

\ j ij/max 

diagonal term; the minimal expected value of x. is then 


1/2 . 

(b..) 

11 


or 


E (y.) 

^ '^i> “ibi- 

11 


(4-20) 
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Using the definition ,4-6 we can calculate the second expected value required in 


Equation (4-4), since 



is normaly distributed. 



1/2 1/2 

The maximum terms of the v.. • b.. b,. are when i = i so that 

1 ] 11 ]] 


E 


/v., y.\ = V.. E (y.) 

I ij jj, 11 1 


(4-21) 


'max 


The condition number is then derived using Equations (4-4), (4-20), and (4-21) 
as 


C. = b.. • V.. 
i 11 11 


(4-22) 


The equation noted above indicates that the loss of significant numbers of digits 
in the solution can be computed for each element in the solution vector using 
the product of the diagonal elements in the original and inverse matrices. The 
error analysis holds only when the solution vector has been iterated such that 
Equation (4-1) Is maintained to the precision of the original matrices. 
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A method of ensuring the accuracy indicated by C. is to iterate the solution until 
the original Equation (4-1) is satisfied. Let 


X. = X . + Ax. 
1 oi 1 


(4-23) 


where x . is the intermediate solution, and x. is the complete solution which 

Ol 1 

satisfies Equation (6-1) . The vector Ay can be computed such that 


Ay = y - Bx^ 


(4-24) 


and x^ is then recomputed so that 


X = X + VAX'- 
0 0 


(4-25) 


After enough iterations x^ approaches the solution vector x . Convergence 
occurs when the ratio 




(4-26) 


is less than 10 , where d is the number of digits accuracy in the matri- 

ces B and y . 
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SECTION 5 - RESULTS 


Extensive tests were made of the generalized Gaussian elimination algorithm. 
The primary concerns were to determine optimum ratios of core requirements 
and computing time. There are many factors that can affect the computing 
time which are functions of the machine that is being used. However, there 
are some variables that can be extrapolated from one machine to another which 
are functions of the algoritlim. These variables are discussed in this section. 

The primary test was to invert a matrix dimensioned at (1500 x 1500) using 
five different partitioning levels. The test was controlled in that the solution 
vector was known a priori. The test showed that numerical stability is not 
affected by the number of partitioning levels. The results also showed that the 
central processing tmit (CPU) time was only moderately affected by the niunber 
of levels of partitioning, whereas the input/ output operating (I/O) time was lin- 
early affected. Figure 5-1 illustrates the results as obtained on the IBM 
S/360-95 compute!’. Figure 5-2 illustrates the change in the l/O slope when 
different sized matrices are Inverted. 

In Figure 5-3, the l/O time is plotted versus the matrix dimension. Core is 
given in units of 1024, K bytes. The plots were obtained using the various par- 
titioning levels. The plot illustrates the problems of inverting very large mat- 
rices, using small amounts of computer core. The I/O time in the plots is 
proportional to the dimension to the fourth power. 

Figure 5-4 illustrates the CPU time required for inverting matrices of various 
dimensions. In this test the CPU time is porportional to the matrix dimension 
to the third power. Since the CPU time is only moderately affected by tiie num- 
ber of partitioned segments (Figure 5-1), the CPU time can be easily established 
for the machine being used. 



RUNNING TIME (MINUTES) 




NUMBER OF PARTITIONED SEGMENTS 


THE 15-SEGMENT CASE USED A TAPE ON ONE OF THE SCRATCH UNITS WITH A LARGER BUFFER SIZE THAN 
THE DISK UNIT CASES (i.e.. 32K VERSUS 7K). 


Figure 5-1. 


Running Time Versus Number of Partitioned Segments for a 
Matrix of Dimension 1500 X 1500 (IBM S/360-95) 



I/O TIME (MINUTES) 



Figure 5-2. I/O Running Time Versus Number of Partitioned Segments 
for a Matrix of Various Dimensions (IBM S/360-95) 


TIME (MINUTES) 





MATRIX DIMENSION 

♦THIS CASE USED A LARGER BUFFER ON ONE OF THE SCRATCH UNITS 
(r.e., 32K INSTEAD OF 7.4K). 


Figure 5-3. l/O Time Required for Inverting Matrices 
Up to 1500 X 1500 


200K 


320K’ 

540K 


1000K 
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Comparisons of running were made for computing the complete inverse and for 
stopping after the forward-elimination procedure and obtaining a solution. 

Table 5-1 gives results of a regression analysis of the resultant CPU and l/O 
curves. The parameters in the table are for inverting a matrix using lOOK core. 
CPU time was fit to the curve 


^CPU ^2 ^ ^ ^2 ^ ’'"^3 ^ 


(5-1) 


and the l/O time 


2 4 

t / =a N+a N +a, N 
I/O 1 2 4 


(5-2) 


The CPU comparisons show that the a^ coefficient is reduced by approximately 
one-third when the solution is obtained without obtaining the complete inverse. 
Likewise, the a^ coefficient is reduced by approximately one-fourth. 

Table 5-1. Coefficients for Computing Running Time 
on the IBM S/360-95 


Coefficient 

a 

1 

a a 

2 3 

a 

4 

CPU (complete inverse) 

2.02 X lO"^ 

1.51 X10~'^ 1.58 XlO'® 


CPU (solution vector only) 

4.76 X 10“^ 

2.27X10"'^ 0.52X10"® 

- 

I/O (complete inverse) 

-1.81 XlO'^ 

3.85X10“ 

1.86 X lO"^^ 

l/O (solution only) 

0. 53 X lO'^ 

3.36 X lO""^ - 

0. 47 X lO"^^ 
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APPENDIX 


Presented in this appendix is a FORTRAN listing of the PINV package, PINV 
is the driver routine for the generalized Gaussian elimination algorithm. 



n o o o o 


SUBROUTINE PINVt A,Y ,X,B tLfM, 10»NN»N) 

ifFffu^e'^rfnTxrfir ~ "" ; _ 

DIMENSION A(NN) tBI l )f 10(1) ft (1) »M{ n 

C ' ■ ■ . ^ - 

C THIS IS THE MAIN DRIVE SUBROUTINE USING RECURSIVE PARTITIONING" 

'C "AS A'GENERALIZEO GAUSSIAN ELEMINATION'TO* INVERT' A-LEAST SQUARES 

C MATRIX. 

nc"" THE SOLUTION FOR X IS OBTAINED IN THE LINEAR SYSTEM AX*Y . 

-C" ALSO THE INVERSE MATRIX IS OBTAINED AS A COMPLETE SQUARE MATRIX.' 

c~- - — - 

C THE PROGRAM EXPECTS THE MATRIX A TO BE AVALABLE ON UNIT AO. EACH~ 

C' RECORD CONTAINS 1 ROW OF AfSTARTING WITH THE DIAGIQNAL ELEMENT.' " 

C 

C the INVERSE MATRIX IS WRITTEN ON UNIT Al IN A SQUARE FORMAT WHERr 

C" EACH ROW IS CONTAINED ON I RECORD. _ 

'c 

C ■ SCRATCH UNITS ARE AO »4 1» A2 t A3 t 4A (II - lA ) 

Q • ^ .. 

C A ARRAY IS USED TO STORE PARTITIONED SEGMENTS OF THE MATRIX A 

C NN IS THE AMOUNT OF STORAGE AVALABLE ' ' " 

C ' B ARRAY IS USED AS A TEMPORARY STORAGE t ID ARRAY ALSO 
C' — X - SOLUTION VECTOR , Y - RIGHT HAND SIDE ' ‘ " 

C N - NUMBER OF ROWS IN MATRIX A 

C ' ' . ■ ■ 

c 

C INITALIZE AND CLEAR ARRAYS - - . 

C 

DO 1 1 = 1, N 
B( I )=0. 

1 X( I )=Y( I ) 

11=40 

12=41 . 

13=43 

14=44 , 

COMPUTE NUMBER OF PARTI TIUNEO SEGMENTS NEEDED (IP), 

NUMBER QF COLUMS IN FIRST ROW (M) , AND NUMBER OF ROWS (L) 

IN EACH SEGMENT 

M( I )=N 
C0=NN*2 
DO 2 1=1, N 
MI=M( I) *(M(I )+l)/2 
IFtMI.LE.i,N)m TO 21 


PAGP to 

POOR QUALITY 
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B0=MM)*2+1 ^ _ 

_ " B0=B0/2 L _ "1~’ 

■“L(I)=80 -SaRT(Bd*B'0-C0) ‘ ' ' ‘ " 

2 M(I + 1} = M( n-Ld ) 

1 21 IP=I _ ■_ 

PRINT 22. IP ■ ■ ■■ ■ 

22 F0RMATI1H0.35HTHE MATRIJ( WILL BE PARTITIONED INT0.I3.1X. 

• 6HLEVELS) 

L(lP» = Hnpr ■” ■ 

CHECK NUMBER OF PARTITIONED SEGMENTS » 1» SKIP PARTITIONING" 
IF(IP.EQ.l) GO TO 4 " J" " ■ -~ 

GAUSSIAN ELIMINATION USING PARTITIONEO^SUBSETS 
IMl=IP-l 

' DO 3 I=l,IMl ' ' ■ ' ' ' ’ y ’ " ' 

READ AND STORE A11.A12 PARTITIONED MATRICES 

LL=L(I ) 

MM=M( I ) 

CALL RDA1A2(A,MM,LL,NN,II) 

INVERT All 

CALL SPGEKLL.MM.A.e.NN) 

COMPUTE A22=A22- A21*A11*A12 AND A12*A11*A12 

CALL FWDELMTA.B.B.LL.MM.NN, IU,I1,I2)' 

COMPUTE X1=A11*Y1 , AND X2=X2-A2l*Xl 
NMM=N-MM+1 

CALL FWDSLVIA.B.XINMMI.LL.MM.NN) 

WRITE DECOMPOSED MATRICES ON UNIT A3 IAII.A12) 

CALL WTAIA2(A,MM,LL’,NN,I3I 
C 

3 CONT INUE 
C 
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C READ UST PARTUIOMED SEGMENT 

z — ^ ^ ^ ^ ^ ^ ^ ^ ’ " 

4 HK*M( IP ) 

' ’■'eALt' RUAlA2(A»MM*HH.NNi IX» ■ 

’ INVERT LASr PARTI TIGNEO SEGMENT - • 

^ CALL SPGEI (MM,MM,A,B*NN) 

COMPUTE X 2 »'A 22 *Y 2 " 1 V 7 ! 

■ nmm • n - mm ♦ I • - 

CALL FtaOSLVU»B,XINPM»,MM,MM*NN) 

WRITE A 22 ON UNIT II . I 

■ Il»I2 ■ 

12 = 12+1 * ^ -r 

' IPn2.GT.+2)12«4l ^ ^ ^ 

CALL WRTB 22 (A,B,Mf 1 ,NN.l 1 ) 

IF THERE IS ONLY I PARTITIONED SEGMENT NOMORE PROCESSING" 
IS REaUIRfcU. 

IFIIP.EG.DRETURN 

OBTAIN SOLUTION THROUGH BACK SUBSTITUTION 
I = IP 

DO 5 J-lt IMl 
l»I-l 

READ M 1 .A 12 

MM=M(I ) 

LL=L(l) 

CALL RpAlA 2 IA,MM,LL.NNtI 3 ) 

SOLVE Xl=Xl-Al2*X2 
NMM=N-MH+I 

CALL BCKSLVU.X (NNM)*XfNMM) tLLtMM.NN) 

COMPUTE nil=All*M2^B22*A2l 


A -4 


B12=-A12*822 . 

.1 BCKSUB ( A»B ,B »LL,MM,NN»1dVI It I2f.K‘) 

5~C0NTINUE _ ' ' ■_ 

RETURN ■■ ■ ■ ■ 

ENO 



SUBROUriNt BCKSLV(A,Xl« X2fLfW ,N»NN) 

'blMENSION A{N>4) 

_«E AL *8 X I ( L » t X2 ( Mj _ ^ 

;^SpL VE S X 1^ ^X 1 A 12*^ 

jj=0 ■ ■ _ 

z!^*\ZZ . 

DO 1 I-ltL ■ * ■ ■ 

^ IJ=IJ+LPlHr‘J~‘ 

DO 1 J=LP1»M 
IJ= IJ+1 

l ' Xl( I )=X1( 1 )-A( rj)+X2( J) 

■"'RETURN 
END 
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I . I • ' I ; ■ . I . ' ■ 


SUBROUr INE BCKSUB ( A »8 > 0 > LtHf NN> I D« II , 1 2t 1») 

"TH I S' SUBROUT INE ENL ARGE'NlT'f Hf^'P ARYff fOTTEir^IB'Mfe^^^ 

"INVERSE MATRIX THROUGH THE FORMULAS ’ ' _ 

Bll= All ♦ A124>b22*A2L ' " 

012= -A12*B22 


DIMENSION AINNI ,BIL )»D(P» ' - - . 

iNTEGER*2 lorn ' ^ _ZZZZriI7TI 

A' ARRAY CONTAINS All ANC Ai2"' " 

' ”* 'B ARRAY IS A TEMPORARY ARRAY* FOR '&21 MATRIX 

0 ARRAY IS A TEMPORARY ARRAY FOR "B22 KATRiX* ' 

' ■' FOR EFFICIENT CORE USAGE B AND D SHOULD ^ “ 

CLEAR ~~~ 

DO I i=itL* ■ 7 ziTiz^~ZTzr” 

1B(I)=0. - ...... - _ - .. 

MM=M-L T i ” T* 

] _1 ' 

7" ; initalize Index ARRAY v iNDExti, jr vt.ge.i ~ 

ID(1)=0 ■ • ~ 

DO 2 I = 2fHM ' " ' 

' 2 iD(i)=io(i-n+M-m - - 

READ B22 FROM UNIT II 
COMPUTE Bll= fill + A12*B22*A21 
B12= -A12*B22 

“ 00 5 K=LPi,M 

CALL READ6(D(LPl)»MM,in 
C 

00 3 1=1, L 
DO 3 J=LP1,M 
IJ=ID(I)+J 

3 B( I) = BM )-A( IJ)»0(J) 

C 

DO 4 1=1, L 
IK=II)(I)+K 
00 4 J=I,L 
ij=iD{n*j 
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4 A(IJ)=ACIJ}-AnK)*B{J ) 

C AOr WR ITE"6 CBTUt Rl 

5 CONTINUE ■ ' 

“C - 

REWIND 14 ' ^ 

REWIND ii^" 

"C READ B2 1 A'N0~ST0RH“1N“C0RE~ 


~ DO 6 K=LPIVM 

CALL' READ6(B'»LVr4T 

*'“00 6 i=r,i 

“ lK=I0m+K 
~£"A(IK)=Bm 


3 RITE B 1 1 , B 1 2;;;oT) ocrTPur~aNT T~n r 

“"'DO 10 I=1,L ■ 

00 9 J = 1,M 

IF(J-I) 7,8,8 ~ 

“”7 lJ=ID(Jj+l " 

GO TO 9 

“■8~I J=IO{ I )+J ^ “ 

~9 D( j)=A( ij) TTH” 

CALL WRITES fO,M,T2} ' 

10 CONTINUE J '"7’ 

WRITE 821,822 ON UNIT 12 "-'" 
DO 12 K=LP1,M 

CALL READ6(D(LPI),MP,I1) ■ 

DO 11 1=1, L ' ; 

IK=1DII)+K 

11 D( I) = A( IK) 

'CALL WRITES (D,M,I2) 

"12' CONTINUE ■ ‘ ~ 

■ SWITCH UNITS AND REWIND 

REWIND T1 . . - - 

REWIND 12 
11=12 
J2=I2+1 
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IF< I2.GT.'#2H2=A1 
RETURN _ 

END 
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SUBROUT INE FWJELM(A t B tO t L tM ,NN t IO» ll 1 1 2) 

'■ c ■ ■ - ^ 

"c THIS SUBROUTINE compute S' 22 = A22 - A2l*All*Al2 

C ' J'A12= ' Ail^AlZ 1" _■ 

c _ _ ' _ ' ' " _ J1_I_ _ \ I— 

„ DIMENSION A(NN) ,B'{L ),D( M)_ __ _ _ _ _ 

~' INTEGEK*2 10(1) 

^A' A R R A Y“' CON T A IN S"aTi" AN D A12 ~ I 

B ARRAY IS A BUFFER FOR TEMPORARY 'STORAGE"OF" Al 2=An*Al 2 

U ARRAY IS A BUFFER 'FOR TEMPORARY STORAGE OF "A22_ ' " 

for EFFICENT use of core B( 1) 'and _ 0( n“'SHJJOXtr EQUI’YAXENT 
ID IS AN INDEX ARRAY ' INDEX ( I , J )'= 'I OCn + j" WHERE" GE . T ' 

INITALIZE B ARRAY 

DO 1 1 = 1, L _ 

I B( I )=0. 


INI TALI 2E J^D ARRAYj^^^^^^^\ ' 

■ I D (' 1 ) = 0 “ 7 _" 

MM=M-L ~ 

' DO 2 1 = 2, MM Z_1 

2 ID( I )=I D( I- ly^-M-t +X ■ ■ _ ^ ' 

LP1=L+1 

■LMi=L-i ■ _ ; ' 

READ A22 FROM UNIT II- 

WRITE A22=A22 - A2 l*Al 1*A12 ON UNIT 12 
STORE A12= AU*A12 IN CORE 


DO 7 K=LP1,M ' 

CALL READ6 (B(K),MH,Il) 

DO 3 1=1, LMl 

IJ=TD(I)-H 
lK=IU(l)-t-K 

B( I) = Bn)+A(IJ)*A(IK) 

IP 1=1+1 

00 3 J=1P1,L 

IJ=IDm+J 

JK=ID(J)+K 

B(I )=B( I)+A( IJ)+A( JK) 
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3 B(J)sB(J)4^A(IJ)*A(IK) 

rj=TD(r)vc — 

IK*10(L)+K 

' B(i )=B( u+At nr^Ad K) 

c 

0GT5 I*KVK' "" 

■ DO 5 J=1,L 

IJ=ID( J)+l ■ 

"5 D| I ) = D { n - A ( IJ) *B ( J )~ 

c 

'DO 6 J= l,L ■■ 

JK=IO( J)+K 

A( JK)=B( J) 

■ ■ “ 6 B{ J} = 0. ~ 

CALL WRITE61D(KJ;HH,I2) ~ 

7 

'C, ' ■ ■ ■ 

REWIND II • 

■' ' REWIND 12 

M=I2 

— * 12 = 12+1 • ’ 

I F ( I 2.5T.42 1 TZ=4I " 

RETURN 

•■■■ END '• 
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SUBROUTINE FHOSLVt A *3 ,X ,L t M,NN) 

irrfTETgsiowAiNni 

REALMS X(M) ,8(L) 


SOLVE XUAll^Xl 

■X2=X'2-A2l*Xl ~ _ 

LP1=C+1 ” 

_ ij=o 

DO A 1=1 a " 

- ij=rjvi 

B( I )=BnT+A(rJl»XU )" 

ipi=i+i 

IF( IPi.GT.DGO TO '2 “ 

c 

DO 1 J=IPliL " ■ 

■ ‘ IJ = IJ + 1 

B( J)=B( J)+A( IJ )*X( 1) 

1 ‘CONTINUE 

~(T ■■■ 

2 IF( LP1.GT.M)G0 TO 4 


DO 3 J=LPltM ■ 

IJ=IJ+1 ■ 

X( J)=X( J)-Al IJ)*Xl I) 

3 CONTINUE 
C 

■ ‘ - x( i)=Bm 

■ B( I )=0. 

4 CONTINUE 

RETURN 

-■ " - end 
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SUBROUTINE R0A1A2CA tM,L tNNt 1 1) 


READ All AND_A12_INTQ ARRAY 3 t READ UNIT II 
DIMENSION A(NN) 


IF A IS READ ON UNIT A3- THEN BACKSPACE THf^OUGH L ROWS 

’MF( I1.NE.43JG0 TO 2 

DO 1 I = 1,L ■ ■■ 

1 BACKSPACE II 


READ All,A12 

2 Hl=l . 

MM=M ■ ■ 

DO 3 1 = 1, L 

CALL READ6{A(HI);MH',TIT 

Ml=Ml + MM ■ ~ 

3 ■■■■ '• 

■c" 

IF( I l.NE.A3)RETURN 

■ — DO A 1= 1,L • ■ 

■ “'A- backspace n --r •• “ 

• " RETURN ■ ■ 

■ • END 
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SUBROUTINE SPGE I(N, M,A,B,NN) 


~z 

C ~ " M ATRI X INVERSIO N US ING"GAPS S1AN~ ELI MATT' O N 

C A ARRAY IS UPPER TRIANGLE OF A POSITIVE OEF I NI TE t SYHETR IC MATRIX 
""C ARRAY I S’ FOR ‘■TEMPORARY“5TORAOE“ OF~BT2~PARTmONEO ’T^OW 

DIMENSION ATNN),B(N) — 


C33INDEX FUNCTION. INDEXC I7J > ~^~I NDn V * J 
IND( I) = (M*2-I)*n-r)/2 — 


T M - NUMBER OF ELEMENTS IN“THE~FIRST KOW^F A 'HATRIX 

C N’- NUMBER OF ROWS IN A MATRIX’ TO BE INVERTED 

"C'~J" NN- NUMBER OF ■ELEMENTS’"IN”A MATR IX t~~ NN= ~ N*H ~-N* ( N-DT2 

C FORWARD ELIMINATIONt All=’'l/A’l’l 

C ■ A22=A22 - A2l*All*AT2' ' 

*^Tn~ ' ' " A’12=AT1*A12 

DO ’2 I = 1 , N H I 

' iPi-n-1 ’ ■ 

■ lNOI«INDm - ”*• 

1 1= INDI + I 

A(II)=1./A( II) “ 

C ■ ' 

00 2 J=IP1,N • ■ ■ 

INUJ=IN0(J) ’" ■ ' 

IJ=1N0I+J • • ■ ■ 


DO I K= J,N 

JK= INDJ + K" • — — 

IK=INOI+K - 

j__- j JK )-A ( I J)'*A ( I D *A riKT 

"C 

2 AnJ) = A(II)<'A(IJ) 

■'c ■■ 

■ ■ tI=rND{N)+N 

A(in = i./A(ii) 

■ C, ■ '• 

C BACK SUBSTITUTION, B22=A22 



ORIGINAL PAGE IS 
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B12» -A12|^22 
fall="All - ai2*A2r 


"CLEAR B ARRAY* 

~ * DO 

>*o. 

I'^N * ■■ 

00 *8 I r=2 1 N 

” " Jl-I * 

-1 = 1-1—'— - 
iNor^iMOTri 


CO^MPUT E B 1 2 A 12 * B 22~~ Wff~K7ge7'J 

■■""""DO 5 J=J1,N 

INDJ=IN0( j) * 

DO* 5 K=J,N ' 

■IK=IN01+K ■ 

JK= INuJ + K 

5 B{ J l=B( J»-A( IK)*A( JK) ■ 


COMPUTE B12 FOR "K.LT.J 

' IF( J1.EQ.N)G0 TO 7 
J2=J1+1 


■ DO 6 J=J2,N 

Kl=J-l 

JK= INUI+J ■ 

K2=I1 

00 6 K=JI,K1 
K2=K2-l 

IK= INOI+K ' ■ ”* 

JK=JK+K2 ■ 

■ 6 B{ J)=(3 U) - A( IK)*AUK) 

7 11= IMDI + I 

COMPUTE Bll = Bll - A12*B2l 

DO 8 J=J1,N 
IJ=IN0I+J 

A(I I) = A( II) - B( J)«A( IJ) 
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A(IJ)=B( J) 
1 BU)=0. 

RETURN ~ 

END 






'/AgjS) to 
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SUBROUTINE WTAl A2 ( A ,M , L tNN, I 1 



C __WRITE All,A12 ON SCRATCH' UNIT 1 1 T~'RD¥S~VH-C0LLUHNS~ 

MM=M 

Mi=i ■■ ‘ 

no 1 1=1. L .. .. .. ~ 

CALL WRITE6(A(MU ,MM,in 

.... m1=MI+MM ' ~ ' 

1 NM=MM-i '■ ■ ■ “ ~ ■ ■ " 

RETURN " 

■ END ■■ 
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